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ABSTRACT 


The  most  used  methods  for  solving  algebraic  eigen-value  problems 
are  discussed  in  detail  in  this  report.   These  methods  are  Jacobi's, 
Householder's  and  QR-algorithms.  Attempts  are  made  to  adapt  or  even 
modify  the  algorithm  to  take  advantage  of  parallel  computations. 

A.  Jacobi's  Method 

It  proved  to  be  one  of  the  most  effective  methods  on  a  parallel 
computer.  However,  it  is  limited  to  symmetric  matrices  with  dominant 
principal  diagonals.   Evaluation  of  the  eigen -values  and  eigen-vector 
is  rather  straightforward  once  the  method  converges. 

B.  Householder's  Method 

This  method  reduces  the  symmetric  matrix  to  the  tridiagonal 
form  by  means  of  elementary  Hermitian  orthogonal  matrices.   Once  the 
matrix  is  in  the  tridiagonal  form,  evaluation  of  the  eigen-values  and 
eigen-vectors  is  performed  by  computational  techniques  that  make  almost 
full  use  of  parallelism. 

C.  QR -Algorithm 

The  QR-algorithm  is  used  for  finding  the  eigen-values  of 
unsymmetric  matrices.   For  best  results,  the  matrix  is  first  reduced 
to  an  upper  Hessenberg  form,  using  Householder's  method.   Convergence 
is  accelerated  by  performing  a  single  origin  shift;  however,  when  the 
matrix  has  some  complex  conjugate  eigen-values,  a  double  origin  shift 
is  used  in  order  to  avoid  complex  conjugate  shifts. 

Flow  charts  and  Tranquil  language  programs  for  all  algorithms 
are  included  in  the  report. 
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1.      Jacobi's  Method 

1.1  Introduction 

There  are  many  methods  for  finding  eigenvalues  and  eigenvectors 
of  a  symmetric  matrix.   Two  of  the  most  popular  methods  are  Jacobi's 
method  and  Householder -Givens  method.   In  order  to  apply  them  on  a 
parallel  computer  effectively,  we  have  made  several  modifications  to 
both  of  the  methods.   In  this  section,  we  will  discuss  Jacobi's  method. 

1.2  Modification  to  the  Classical  Jacobi's  Method 

Instead  of  using  the  Classical  Jacobi's  Method,  i.e.,  going 
through  all  the  trouble  to  search  for  a  pair  (since  it  is  symmetric) 
of  largest  off  diagonal  element  and  to  eliminate  them,  we  modified 
the  method  by  eliminating  all  the  off-diagonal  elements  of  each  2x2 
submatrices  along  the  diagonal  through  orthogonal  transformation. 
At  each  iteration,  this  modification  will  remove  n  elements  of  a 
n  x  n  matrix„   We  eliminate  the  elements  by  determining  an  orthogonal 
matrix  cp  such  that  cpA1  cp  =  A    where  A   ~  having  zeros  in  the  approp- 
riate locations.   To  eliminate  those  elements,  we  choose 

cp  =  Diag.  (Tx,  T2,  ...,  T^2) 

[cos  a       sin  a    ~| 
K  =  1,    2,    ,    n/2 
sin  a,    -cos  a                             '      '  '     ' 

t  t  t 

So  T   =  T    ,   Cp  =  cp     and  TT     =  I. 

If  matrix  A     is   partitioned  into  2x2  blocks   as   follows: 


-  1  - 


A1  = 


All      A12      A13       *  '  '   Al  n/2 
, j , r _ 

A12      A22      A23       '  *  *  !  A2  n/2 
r- r p -----r_ - 

At      !  A*      !  A  A 

H13    |   23    ;   33    j   •  *  *  i   3  n/2 

1- 1 1 1 

•  I  •        I  •        I  la 

•  I        •      I        •      I  I        • 

•  •      I        •!•••!        • 

•  I       •     I       •     I  I       • 

f. , , , __ 

Al  n/2  !  A2  n/2  j  A3  n/2  j   '  '  ■  jAn/2  n/2 


Then  A1+  =  CpA1  cp  will  be 


i+1 


T1A11T1 

TAT 
1   12    2 

TlAln/2Tn/2 

T  At  T* 
2   12   1 

T    A      T 
2   22   2 

T2A2n/2Tn/2 

Tn/2Aln/2Tl 

T        A*        T* 
n/2A2n/212 

• 
• 

Tn/2An/2n/2Tn/2 

Considering  the  diagonal  siibmat  rices  A,  ,      =  T.  A, .  T, 
to  to  kk  k  kk  k 


cos  OL       sin  a, 
sin  a,    -cos  OL 


a2k-l,2k-l     a2k-l,2k 
J*2k-l,2k         a2k,2k    J 


cos  Oi       sin  cc 
sin  <X    -cos  OL 


b2k-l,2k-l     b2k-l,2k 
_b2k-l,2k         b2k,2k    J 


where     ^_1>2^_±  u  *2k-l,2k-l  cos~< 

+  a2k,2k  Sin2  \> 


k       "2k-l,2k  Sin     °k 
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b2k,2k  =  a2k-l,2k-l  Sin2  ak  "  a2k-l,2k,  ^  \ 
+  a2k,2k  C°s2  \ 


and  b2k-l,2k  "  2  (a2k-l,2k-l  '  a2k,2k] 


sin2  ak   -  a^      cos2  ^ 


Now  to  eliminate  the  off -diagonal  terms,  i.e.,  b       =  0 
tan  ak  =  (a2k_1;2k)/2  (a2k-l,2k-l  "  a2k,2k) 


1  x       "I 
a     -  —   tan 


2  a 


2k-l,2k 


^a2k-l,2k-l  "  a2k,2k> 


If  we  define  that 


z   _  I 

Zk  ~  2  \ 


1/([Ua22k-l,2k/(a2k,2k  "  a2k-l,2k-l)2]  +1^ 


Then  the  rotation  angles  will  be 


sin  <X    = 
k 


—  +  Z.       ;    and   cos  a, 

2k7  1 


\2 


i-z. 


and 


^k-l^k-l     "     2    (a2k-l,2k-l  +   a2k,2k)   - 
2k      ,2k 


\ 


2  /a2k-k,2k-l       a2k,2kx2 

a2k-l,2k  +    [  2  ; 


we  also  notice  that 


i+l> 


tr    (A1)    =  tr    (A1+±) 


and  since  the  norm  is  invariant  under  an  orthogonal  transformation  of 
the  matrix.   Therefore, 

N2  (A1+1)  =  H2  (A1) 


-  3  - 


Since 


N2  (A1)  =  tr  (A1  A1) 


where  A  is  the  absolute  of  A.  Also 


2  2        2  2  2 

therefore, b2k-l,2k-l  +  b2k,2k  =  a2k-l,2k-l  +  a2k,2k  +  2  a2k-l,2k 


Hence  by  one  orthogonal  transformation,  the  sum  of  squares  of 
the  diagonal  coefficients  increases  by  the  sum  of  squares  of  the 
(n)  vanished  terms. 

Now,  after  one  orthogonal  transformation,  the  matrix  will 
look  as  follows: 


A 


i+1 


\i°    !  bi3  \h\-  ■ — |vv 

0       b22      "      b23   b2U -!b2,n-lb2,n 

J-                       1                        1 

h13h23    j    h33    o    j ; 

\kh2k      1      °       \k   I"     "      "      -' " 

.L                       1                        i. 

J_                       1                        J_ 

_l_       _       _       _l_       -       -       -lh                        0 

1                                               j    n-l,n-l 

-■-       -       -       -i-       -       -       -0               h 

1                        '                        '                 n,n 

Therefore,  "by  shifting  the  second  row  to  the  "bottom  and 
the  second  column  to  the  far  right,  we  bring  to  the  diagonal  new 
submatrices  with  non-zero  off-diagonal  elements,  and  hence  the 
matrix  is  ready  for  the  second  transformation  .   This  shifting  is 
repeated  for  n  -  1  transformations.   Then,  we  will  shift 


-  h   - 


the  first  row  to  the  bottom  and  first  column  to  the  far  right 
repeat  the  above  mentioned  shuffling  process  for  m  transformation 
(the  diagonal  terms  are  sufficiently  dominating).  Therefore,  the 
diagonal  elements  will  be  equal  to  the  eigen-values,  since  the  eigen- 
values are  not  changed  by  orthorgonal  transformations.  For  comparing 
number  of  transformations  required  by  this  revised  method  with  the 
classical  Jacobi's  method,  see  Appendix  F,  first  Q.P.R.  The  eigen- 
vectors will  be  the  column  vectors  of  B,  where 


B  =  <Plt  *2      •••  C 


because  if  we  let  the  final  matrix  be  F,  then 


t       t  t       „ 

cpm  . . .   cp3  cp2  cp1  A  cpx     cp2      . . .   9m     =  F 


but  AX  =     XX 

Since  F  is  almost  diagonal  with  X's  being  its  diagonal  elements, 

. * .  AX  =  FX  where  X  is  the  matrix  with  eigen-vectors  as 

its  columns 

t        t  t       _        t       t  t 

Hence  A  cpx     cp2      •  •  •   <f>m     =  F  cp1     cp2      ...  <Pm 

and  B  =  cp.,   cp^   ...  CD    is  the  matrix  with  the  corresponding 

Tl  ^2      Tm 

eigen-vectors. 

Thus,  in  practice  one  can  follow  the  following  four  steps: 

(i)   Calculate  orthogonal  transformation  matrix  cp  by 
finding  n  rotation  angles  and  construct  2X2  submatrices  on  the 
diagonal  and  zeroes  everywhere  else. 


(ii)   Pre  and  post-multiply  matrix  A  at  i-th  step  by 
.  transformation  n 
;et  B    by  B  cp  where  B  ,  initial  B,  is  identity. 


orthogonal  transformation  matrix  cp  to  get  A   .   At  the  same  time, 


-  5  - 


(iii)   Testing  whether  the  ratio  between  the  sum  of  the  squares 
of  the  diagonal  elements  and  the  sum  of  the  squares  of  the  off -diagonal 
elements  is  greater  than  (l/£),  where  £  is  a  given  small  number  (.001  - 
.01).   If  greater,  print  out  the  eigen-values  and  the  eigen-vectors  from 
the  present  matrices  A  and  B  respectively,  otherwise  go  to  next  step. 

(iv)   Shuffle  the  second  row  and  second  column  to  start  another 
transformation.  After  (n-l)  such  shuffling,  shuffle  the  first  row  and 
first  column  and  start  all  over  again.   If  after  (n-l)  shuffling  of  the 
first  row  and  column,  the  diagonal  elements  are  still  not  sufficiently 
dominant,  the  method  fails. 

Jacobi ' s  method  has  proved  to  be  suitable  for  diagonally- 
dominant  matrices.   Therefore,  although  Householder's  method  works  for 
any  symmetric  matrix,  due  to  the  simplicity  of  Jacobi' s  method,  it  will 
be  used  for  finding  eigen-values  and  eigen-vectors  for  diagonally-dominant 
matrices. 

1.3  High-Level  Language  Code  and  Flow  Chart 

For  flow  chart  and  high-level  language  code  see  ILLIAC  IV 
Document  Number  165. 

2.   Householder's  Method 
2.1  Introduction 

This  method  is  one  of  the  most  effective  methods  that  exist 
for  reducing  a  symmetric  matrix  to  tri-diagonal  form  using  elementary 
Hermitian  orthogonal  matrices.   For  a  matrix  of  order  n,  there  are 

(n-2)  steps  in  this  reduction,  in  the  r   of  which  [n-(r+l)]  zeroes 

th  th 

are  introduced  in  the  r   row  and  r   column  without  destroying  the 

zeroes  introduced  in  the  previous  steps. 


-  6  - 


2.2  Theoretical  Background 

The  present  problem  is  that,  for  a  given  vector  x  of  order 
n,  it  is  required  to  determine  an  elementary  Hermitian  matrix  P  such 
that  pre-multiplication  by  P  leaves  the  first  component  of  x  and 
eliminates  the  rest,  i.e., 

Px  -  ke  (2.1) 


where 


M 


constant 


e  =  (  1,  0,  0,  ....,  0}   (column  vector) 

Let  P  =  I  -  2wwt  (2.2) 

where  w  is  a  column  vector  of  order  n  and   w  w  =  1.0.   P  is  not 
only  elementary  Hermitian  but  also  orthogonal  for 

P*  P  =   (I  -  2wwV  (I  -  2wwt) 
=  I  -  WwJ  +  Uw(w  w)w 

=  I 


from  Equ.  (2.1),  it  can  be  proved  that  the  Euclidean  norm  of  x  is 

invariant,  i.e. 

02     2     2         2    2 
S  =  x1   +  x2   + xn  =  k  , 

f  t 

2-Norm  of  Px  =  (Ex)   (Ex)  =  (Ex)  (Ex) 

=  xt  P*  F  =  xVpx 
x 

=  x  x  which  is  the  2-norm  of  x, 

Therefore,  k  =  +  S  (2.3) 

hence,  equation  (2.1)  yields, 


x   -  2w  K  =  +  S 

x.  -  2w.  K  =  0   (i  =  2,3, ...,n) 


(2.U) 
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where  K  =  w  x  which  is  a  scalar.   Squaring  and  adding  equations  (2.4), 
the  value  of  the  scalar  K  can  be  obtained 


K 


\I77 


x.  S 


(2.5) 


If  u  =  (x  +-  S,  x  ,    x„,  ....,    x  ),  equations  (2.k)   yield, 


w  = 


2K 


and 


P  =  I 


uu 


2K 


(2.6) 
(2.7) 


For  a  matrix  A  of  order  n,  suppose  that  the  configuration  of  the  matrix 
just  before  the  rth  step  is  as  follows: 


A 


r-1 


1 

0 

c 

r-1 

1 
1 

1 

^Vx 

0 

T 
1 

b 
r- 

_  T 

■l1 

r-1 

1 — 

1 

1 

(2.8) 


where  C  n  is  a  tri-diagonal  matrix  of  size  r,  b  ,  is  a  column  vector 
r-1  r-1 

of  order  (n-r)  and  B  n  is  a  matrix  of  size  (n-r).  The  matrix  P  of 
the  rth  transformation  may  be  expressed  in  the  form, 


"    I      »    0      " 
1 

'     I      '      0                  "1 

I 

.   o    '  Qr   _ 

0      •      I-2v  v  t 
l.                           r  r    -* 

(2.9) 


where  v  is  a  unit  vector  having  (n-r)  components,  w   =  (0  |  v  ). 
Hence,  we  have 


A  =P  A  _  P 
rs  r  r-1  r 


'r-1 


_  _  r  _ 

0  i  C 


r  r-1  r 


(2.10) 
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where  c  =  Q  b  , 
r    r  r-1 

If  we  chosse  v  such  that  all  components  of  c  vanish  except  for  11 
first  component,  then  A  is  a  tri -diagonal  matrix  as  far  as  its  first 
(r-t-l)  rows  and  columns  are  concerned. 

Making  use  of  equations  (2.2)  through  (2.7),  the  method 
can  be  organized  in  a  fast  systematic  manner, 

P  =  I  -  (u  u  t/2K  2) 


where 


where 


u.r  =  0     (i  =  1,  2,  ....,  r) 

u   ^  =  a   J_1  ~  S  (2.11) 

r,r+l    r,r+l    r 

u.   =  a  .    (i  =  r  +  2,  . . . . ,  n) 
lr    ri 


n       1 

S   (   E  a^  .)2 
r   .    .   ri 
i=r+l 


1 


K  =  |"(S  2  -  a   j_.  S  )/2]2 
r      r  -t-   r,r+l  r  ' 


In  equations  (2.11)  the  sign  associated  with  S  is  to  be 

taken  as  that  of  a    , . 
r,r+l 

The  method  discussed  above  involves  quite  a  large  number  of 
multiplications  which  may  be  necessary  in  case  of  reducing  a  non- 
symmetric  matrix  to  an  upper  Hessenberg  form.  However,  for  symmetric 
matrices,  the  above  method  can  be  modified  slightly  in  order  to  take 
advantage  of  the  present  symmetry. 

A=PA,P  (P  t  =  Pj 

r    r  r-1  r  r     r 


using  equation  (2.7),  and  denoting  that, 


-  9  - 


(A     .    u   )/2K  2  =  p        and      (u  t  A     J/2K  2     =     p  t  (2.12) 

r-1     r  '      r            r                     r        r-1''       r                 r  v          ; 

.'.      A     =  A     -    -  p     ut-u     pSu      (u  t   p   )    u  t/2K  2  (2.13) 

r          r-1         r     r            r     r            r       r       r       r  '      r  ' 


again  defining, 


then, 


4 .  =P.  "|u   (u*  p  /2K2)  (2.14) 


r   ^r   2  r    r   r   r 


Ar  -  \_i  -  %  \     -   (qr  \^  (2-15) 


The  matrix  of  interest  is  (Q  B    Q  ),    Equ.  (2.10),  which 

is  of  order  n-r,  since  the  elements  in  the  first  (r-l)  rows  and 

columns  are  the  same  as  the  corresponding  elements  of  A  n.      Therefore, 

r-1 

the  transformation  matrix  P  is  designed  as  to  eliminate  a  .  (i  =  r-f-2, 

. ...,n)  and  to  convert  a    _  to  +  S  ,  where  S  is  defined  before. 

r,r+l    —  r'        r 

Hence  the  vector  p  has  its  first  (r-l)  elements  equal  to  zero,  Equ. 
(2.12),  and  the  rth  element  is  not  needed.   Therefore  there  are 
(n-r)   multiplications  required  for  the  evaluation  of  p  .   Similarly, 
it  can  be  proved  that  there  are  (n-r)  multiplications  involved  in 

the  computation  of  q  .   Finally  the  total  number  of  multiplications 

r  2 

involved  in  the  computations  of  the  elements  of  A  is  [2(n-r)   + 

(n-r)],  but  since  the  matrix  is  symmetric,  the  number  of  multiplica- 

tions  concerning  the  upper  triangle  only  is  (n-r)   +  (n-r)/2,  the 

second  term  may  be  neglected  compared  to  the  first,  .'.   No.  of 

2 
multiplication  for  one  transformation  ~  (n-r)  .   Hence,  the  total 

number  of  multiplications  to  reduce  the  matrix  to  a  tri-diagonal  one 

2  3 
is  essentially  rr  n  which  is  half  of  that  of  Givens'  method. 

2.3  Computing  Eigen-Values  and  Eigen-Vectors 

2.3.1  Eigen-Values 

Once  the  matrix  has  been  reduced  to  the  tri-diagonal  form, 
there  are  so  many  techniques  for  getting  the  eigen-values,  for  example 
(Bisecting  method,  Sturm  sequence  or  Gerschgorin  discs). 
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An  effective  method  for  computing  eigen-values  that  are  most 

suitable  for  parallel  computers  can  be  summarized  as  follows.  The 

largest  and  smallest  eigen-values  can  be  obtained  by  the  power  method. 

Dividing  the  interval  between  \  .   =  \.    and  \         =  \     into  (n-l) 

nun    1      max    n         ' 

divisions  and  get  \.    (j=l,2,....,n).  For  each  \.,   the  determinate 

J  J 

of  the  tri -diagonal  matrix  can  be  obtained,   using  the  method  of  prin- 
cipal minors 


f      (\.)    =  1 
o      y 

fl   ^j)   -  Can  -  V,) 


(2.16) 


fl   <XJ>    =   ^Sk-  V    fk-l(^}    "   Vl,k  WV 


supposing  that  the  value  of  each  determinate  corresponding  to  X .  is 

J 
R.  ( j=l,2, . . . . ,n)  a  plot  can  be  made  between  \.  and  R.,  Figure  a. 


R 


T 


mm 


Every  R=0  indicates  the 
position  of  an  Eigen-value 


max 


Figure  (a) 
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Each  zero  crossing  indicates  a  correct  eigen-value,  supposing  that 

X.    gave  use  to  a  positive  residual  R.,_  while  X.  gave  a  negative 
J+l  J+l        J 

R.,  then,  again  the  interval  between  \.    and  \.        can  be  divided  to 

J  J      J  -^- 

get  the  exact  eigen- value  with  enough  accuracy. 

2.3-2  Eigen-Vectors 

If  x  is  the  eigen-vector  of  the  tridiagonal  matrix  (A   ) 
corresponding  to  an  eigen-value  \,    then  the  corresponding  eigen-vector, 
V ,  of  the  original  matrix  A  is  given  by, 


V  =P1P2  ••■  P„-2X 


To  avoid  storing  the  matrices  P.,  the  vector  v  is  calculated 
from  the  relations, 


P  0  x  =  (x)   0 
n-2     v  7n-2 


Pn_3(x)n-2  =  Wn-3 
v  =  p1(x)2  =  (x\ 


(2.17) 


So,  it  remains  to  find  the  eigen-vectors  of  the  tri-diagonal 

SS' 

in  the  form, 


matrix.   Assume  that  the  tri-diagonal  matrix  (A  „)  can  be  expressed 

n-2 


n-2 


cl  1 

bl  i 

0 

0 

0 

0 

0 

0 

L               i 

\ 

0 

0 

0 

0 

0 

cs 

• 

0 

0 

0 

0 

. 

. 

0 

0 

0 

0 

• 

c. 

1 

b. 

l 

0 

0 

0 

0 

b. 

1 

ci+i 

b.    , 
l+l 

0 

0 

0 

0 

bi+i 

• 

• 

0 

0 

0 

0 

• 

• 

0 

0 

0 

0 

0 

# 
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for  an  eigen-value  X   and  an  eigen-vector  x,   where  x     =   (x    ,  x0,    ..., 
x    ) ,    of  (A     ,,)    the  following  relation  can  be  obtained, 

(C±   -  X)  x±  +  b1  x2  =  0 

Vl  xi-l  +   (Ci  "  X)  xi  +  bi  xi+i  =  °  U=2,3,...,n-1) 
Vl  xn-l  +   K  "  M  xn  =  0  (2.18) 

If  x      =  0,    the  rest   of  the   components   of  the  vector  x   are   zero,    so 
x     will  be   assumed  as   1.      From  the  equations    (2.18)    and   (2.l6)    an 
expression  for  x      (r=2,...,n)   may  be   obtained  in  form, 

xr  =   (-l)r_1  f  (X)/(b1  b2   ...   br)    (r=2,...,n)  (2.19) 


Similarly  all  the  eigen-vectors  of  the  original  matrix  A  may  be 

o 

obtained,  each  corresponding  to  an  eigen-value. 

A  flow  chart  and  a  higher  level  language  program  for  this 
algorithm  to  find  all  eigen-values  are  in  sections  2.k. 
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.4     High  Level  Language  Program 


#BEGIN 

#LABEL 

#ARRAY 


START 


L00P1 : 


#REAL 

#REAL 

#SET 
#READ 


START,   L00P1,   LOOPE,   LP,   LLP,   LI,   LL1,   L2,    SEUP,   REPE, 

FINAL,    FN>L,    EEND; 

#REAL     #FLOAT     #SHORT     #STRAIGHT 

A[  614,64],   MQ[64,64],   TMQ[ 64,64],    FNT[64,64],   FCN[64,64], 

UT[64,l],   U[l,64],    P[l,64],    Q[l,64],    EVALUE[64], 

FNTO[64],    FNTl[64],    FTO[64],    FTl[64],   NEVAL[64]; 

#FLOAT     #SHORT 

S2,    S,    T,    TEM,    GRAG,    INTV,    STRT,    FINI,    RANG,    STP, 

SSTP,    INIA,    FINB; 

#FIXED 

J,    CNT,   NCNT,    EV,    TK,    JK,   NUM,   NUM1,    H,    L,   K,    I,    N,   M; 


HH,    LL,   KK,    II,    JJ,    IJ,    JI, 

A; 


I,   NN; 


J  *-  2; 

HH    : :    [J,J+1,    ,    64]; 

LL    ::    [1,    2,    ....,    J-l] ; 

KK    : :    [J+1,    J+2,    ,    64]; 

S2  <-  #FOR   (H)   #SLM   (HH)#DO  #SUM(A[  J-l,H]t  2) ; 

S  <-  SQRT   (S2); 

T  ^  S2   -  A[J-1,J]    x  S; 

#FOR   (L)   #SIM  (LL)#DO       UT[L,l]   -  0.; 

UT[J,1]   -  A[J-1,J]    -   S; 

#FOR   (K)#SIM  (KK)#DO  UT[K,l]  -  A[j-1,K]; 

U  «-  TRANSPOSE [UT]; 

P  ^  A  x  U; 

P  -  P/T; 

TEM  <-  UT  x  P; 

TEM  *-  TEM/T; 

Q  *-   .5  x  TEM  x  U; 

Q  -  P-Qi 
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MQ  «-  Q  x  UT; 
MQT  «-  TRANSPOSE[MQ]; 
A  -  A  -  MQ  -  MQT; 
J  «-  J+lj 
#IF       J  <  63  #THEN  #GOTO  L00P1; 
NCNT  «-  0; 
EV  «-  1; 

GPAG  «-  EMAX[A]  -  EMIN[A];  (SEE  APPENDIX  A) 
INTV  <-  GRAG/8.  ', 
STRT  <-   EMAX[A] ; 
II  ::  [1,2,... .,64]; 
JJ  ::  [3,h,....,6k]; 
LOOPE:  FINI  «-  STRT  -  INTV; 

RANG  <-  STRT  -  FINI; 
STP  <-  RANG/ Gk.; 

#FOR      (I)  #SIM  (II)  #D0 

#BEGIN   EVALUE[I]  «-  FUNI  +  STR  x  (l-l); 

FNT0[I]  <-  1.; 

FNT1[I]   <-  A[l,l]    -  EVALUE[I]; 
#END; 

#FOR  (K)  #SEQ   (JJ)   #D0 

#FOR  (I)   #SIM   (II)#D0 

#BEGIN 

FNT[K,I]   «-   (A[K,K]-EVALUE[I])    x  FNTl[l] 

A[K-1,K]T2-FNT0[I]; 

FN0[I]   -  FNT1[I]; 

FNT1[I]   *•  FNT[K,I]; 
#END 

NCNT  -  NCNT  +  1; 
TK  «-  1; 
CNT  «-  0; 
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LP :  #IF  TK  >  Gk  #THEN     #GOTO     FINAL     #ELSE      . 

#IF  FNT[6U,TK]    <  0.   #THEN     #G0T0     LI     #ELSE 

#IF  FNT[6U,TK]   >  0.     #THEN     #G0T0     L2     #ELSE     #D0 

#BEGIN 

#PRINT     EVALUE   [TK]; 
EV  «-  EV  +  1; 

#IF  EV  >  63     #THEN     #G0T0     EEND     #ELSE     #D0 
#BEGIN 

TK  <-  TK  +   1; 
#G0T0  LP; 
#END; 
#END; 

LI :  #IF  TK  >  Gh     #THEN     #G0T0     FINAL     #ELSE 

#IF  FNT[64,TK+l]    <  0.   #THEN     #D0 

#BEGIN 

TK  <-  TK+1;  #G0T0  LI; 
#END     #ELSE 
#IF  FNT[6U,TK+l]    >  0.      #THEN     #D0 

#BEGIN 

TK  -  TK+1;  #G0T0  SEUP; 
#END     #ELSE     #D0 
#BEGIN 

#PRINT  EVALUE    [TK+1]; 
EV  -  EV  +  1; 

#IF  EV  >  Gh     #THEN     #G0T0     EEND     #ELSE     #D0 
#BEGIN 

TK  <-  TK+1;  #G0T0  LI; 
#END; 
#END; 

L2:  #IF  TK  >  6^     #THEN     #G0T0     FINAL     #ELSE 

#IF  FNT[64,TK+l]    <  0.      #THEN     #D0 

#BEGIN 

TK  *-  TK+1;  #G0T0     SEUP; 
#END     #ELSE 
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#IF  FNT[6>+,TK+l]   >  0.     #THEN     #D0 

#BEGIN 

TK  ♦-  TK+1;  #G0T0  L2; 
#END     #ELSE     #D0 
#BEGIN 

#PRINT  E VALUE   [TK+l]; 
EV  *-  EV  +  1; 

#IF  EV  >  6k     #THEN     #G0T0     EEND     #ELSE 
#BEGIN 

TK  <-  TK+1;  #G0T0     L2; 
#END; 
#END; 

SEUP  MM  «-  CNT  x  16  +   1; 

NUM1  -  HUM   +  15; 

NN    :  :    [1,2, ,16]; 

MM    ::    [NUM,   WUM+1, ,NUM1] ; 

INIA  «-  E VALUE   [TK-l]; 
FHB  <-  EVALUE   [TK]; 
SSTP  -   (FINE   -   INIA)/l6.    ; 
#F0R    (N)   #SIM   (NN)   #D0 
#F0R   (M)   #SIM   (MM)  #D0 
#BEGIN 

NEVAL    [M]   -  INIA  +  SSTP  x   (N-l); 
FT0[M]   <-  1.    ; 

FT1[M]  *-  A[l,l]    -  NEVAL[M]; 
#END; 
CNT  «-  CNT  +  1; 

#IF  CNT  <  k     #THEN     #G0T0     LP     #ELSE 

REPE  TEM  <-  CNT  x  16 

IJ    : :    [1,2, ,TEM]; 

JI    ::    [3,k, ,TEM]; 

#F0R   (K)  #SEQ   (JI)   #D0 
#F0R   (I)   #SEQ   (IJ)   #D0 
#BEGIN 

FCN[K,I]«-   (A[K,K]-NEVAL[I])    x  FTl[  I]    - 
A[K-l,K]t2  x  FT0[I]; 
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LLP 


LL1 


FINAL 
FNAL: 


FT0[I]   «-  FTl[l]; 

FT1[I]  -  FCN[K,I]; 

#END; 

JK  *-  1 

#IF 

JK  >  Gh     #THEN     #GOTO     FINAL     #ELSE 

#IF 

FCN[6^,JK]    =   0.      #THEN     #D0 

#BEGIN 

#PRINT  NEVAL   [JK]; 

EV  <-  EV  +  1; 

#IF  EV  >  6k     #THEN     #G0T0     EENL     #ELSE     #D0 

#BEGIN 

JK  *-  JK.+  1;  #G0T0  LL1; 

#END; 

#ENL     #ELSE     #D0 

#BEGIN 

JK  -  JK  +  1;  #G0T0     LLP; 

#END; 

#IF 

FCN[64,JK]    =   0     #THEN     #G0T0     LL1  #ELSE     #D0 

#BEGIN 

JK  <-  JK+1;  #G0T0  LLP; 

#END; 

#IF 

CNT   =  0     #THEN     #G0T0     FNAL     #ELSE     #G0T0     REPE; 

#IF 

NCNT  <  8     #THEN     #D0 

#BEGIN 

STRT  «-  FUJI;  #G0T0     L00PE; 

#END     #ELSE     #D0 

#BEGIN 

#PRINT   "SOME  ROOTS  ARE  EITHER  TOO   CLOSE 

TOGETHER  OR  THERE  EXIST  MULTIPLE   ROOTS"; 

#STOP; 

#END; 

EEND :      #END 
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2.5     Flow  Chart 


HOUSEHOLDER  TRIDIAGONALIZATION 


S^   =  S2 


2 
2K      =  T 


H  «-  2 


2        64 
S^-Z  [A((H-1),J)]2 
J=H 


2K2  -  S2-A(H-l,H)x  S 


L«1,2,...,H-1  > 


UT(L)  «-  0. 


UT(H)   <-  A(H-1,H)-S 


K=H+1,H+2,...,N  » 


UT(K)  -  A(H-1,K) 


U  *-  TRANSPOSE (UT) 


[II 


<D 
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p  *-  a  x  u/ac' 


n        td  !        tt        /UTxPN 
Q  «-  P-^-  x  U  x   ( — s-) 


2K 


A  -  A   -  Q  x  UT   -   TRANSPOSE^  x  UT) 


H  «-  H  +  1 


H  <  63 


yes 


<D 


FINDING  EIGEN -VALUES 


no 


FIND  EMAX  AND  EMIN 


DIVIDE  INTO  8  SEGMENTS 
BETWEEN  EMAX  AND  EMIN 


NCNT  «-  0 

EV  <-  1 
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3=1,2,... ,6k 


K=2,3,...,61+ 


NCNT  *-  NCNT  +  1 


th 
DIVIDE  THE  NCNT 

SEGMENT  INTO  6k   INTERVALS 


f  (E.)=l 

f (E  )=(A(1,1)-E.) 

fk(E.)  =  (A(K,K)-E.)xfk_1(E.) 
-A2(K-l,K)xfk_2(Ej.) 


PRINT  THE  EIGENVALUES  FOR  WHICH 
f6[f(E)=0  AND  INCREMENT  EV 

WHEN  EACH  EIGENVALUE  IS  BEING  FOUND 


yes 


STOP 


No 


DIVIDE  EACH  PAIR  OF  E .  and  E    ' s 

FOR  WHICH  ^(E.)  HAVING  DIFFERENT 

SIGH  FROM  fV^E.+l)  INTO  l6  INTERVALE 

AND  DO  FOUR  SUCH  PAIR  AT  A  TIME 


© 
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0* 


No 


<& 


No 


F    (E.)  =  l 
ov    jy 


F1(Bj)=(A(l,l)-Ej) 

FK(EJ)=(A(K,K)-EJ)xFK_1(EJ) 
-A2(K-l,K)xFK_2(EJ) 


PRINT  THE  EIGENVALUES   FOR  WHICH 
f6^(e)=o  AND  INCREMENT  EV  AT 

EACH  EIGENVALUES 


SJ?ECIAL._CASE 
yes 


NCNT  >  8 


yes 


PRINT 
"SOME  ROOTS  ARE 
MULTIPLE  ROOTS" 


AFTER 


"YsTOPj 


H 
H 

O 


o 

H 
1-3 

I 


0^ 

-p- 
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3 •   An  Algorithm  For  Finding  Eigen-values  and  Eigen-vector: i 
For  N on- Symmetric  Matrices 

3 • 1  Introduction 

There  are  two  algorithms  for  dealing  with  unsymmetric 
matrices.   The  first  one  is  (LR)  which  was  developed  by  Rutishauser 
(1953).   The  LR-algorithm  has  proved  to  be  a  powerful  method  for 
finding  the  eigen-values  of  symmetric  band  matrices.   However,  it 
proved  to  be  inferior  for  handling  the  ei gen- value  problem  of  large 
non- symmetric  matrices,  the  most  important  difficulties  are, 

(i)   Triangular  decomposition  may  not  always  be 
possible. 

(ii)   The  amount  of  computation  required  by  the  method 
is  very  great.   The  LR-algorithm  is  essentially 
an  iterative  procedure,  where  the  method  consists 
of  a  series  of  similarity  transformations  on  the 
original  matrix, 

A±  =   Lx  Rx  (3-1) 

where, 

A  is  the  matrix  under  consideration  (size  =  n) . 

L  is  a  lower  triangular  matrix. 
R  is  an  upper  triangular  matrix 

i.e. , 


L  = 


1,1 


'n,l 


where  I.  . 
11 


1,  i  =  1,  n 


-i 


n,n 


r  r 

11      In 


nn 
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by  similarity  transformation, 


hence 
similarly 

hence 

let, 


A2  =   Ll     Ax  Lx 


A2   =   Rl  Ll 


VlBS    ^k\\    ••'    L-1)(A1)(L1L2    ...    Lk) 


(LXL2    ...Lk)   Ak+1=AX    (LlL2    ...    Lk) 


\  =  Ll  L2 


k 


\  =  W-i  •••  Ri 


therefore,   Tfc  UR  =  L  Lg  •••  1^  (i^  \)    \_±    -  •  •  \   ^ 


=  Al  Ll  L2  •  '  *  Lk-1  \-l  ' '  *  R2  Rl 


and  so  forth, 


finally 


T.  Uk  =  <A1> 


(3-2) 


(3-3) 


(3-4) 


and  it  was  shown  that,  if  certain  conditions  are  fulfilled  then,  as 
k  -*  00,  then  (A  )   tends  to  an  upper- triangular  matrix,  where  the 
diagonal  elements  are  the  eigenvalues  in  order  of  their  modulus, 
(the  first  being  the  largest),  i.e., 


(Ax)k  = 


X 


o   \ 


■\ 


n 


(3.5) 


follows: 


The  restrictions  on  the  LR-algorithm  nay  be  summarized  as 

1.   Eigenvalues  must  be  distinct  since  convergence  depends 
upon  the  ratio  (\       /\    )  and  will  be  very  slow  if 
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separation  of  the  eigenvalues  is  poor. 

2.   Triangular  decomposition  of  A  must  exist,  (it  fails 
if  any  of  the  principal  minors  of  A  vanishes). 

3«   The  leading  principal  minors  of  the  modal  matrix, 
(matrix  whose  columns  are  the  eigenvectors),  and 
its  inverse  should  not  vanish. 

h.      Even  if  the  triangular  decomposition  exists,  there  is 
a  large  class  of  matrices  where  triangular  decomposi- 
tion is  numerically  unstable  which  may  lead  to  gross 
errors  in  the  values  of  the  computed  eigenvalues. 

2  3 
5.   The  volume  of  computation  is  very  large,  —  n  multi- 
plications for  each  step  of  iteration,  half  of  them 
in  triangulation  and  the  other  half  in  post-multipli- 
cation. 

Due  to  the  above  mentioned  restriction  on  LR,  QR  proved  to  be 
a  much  more  efficient  algorithm,  and  it  will  be  discussed  here  in  more 
detail. 

A  more  stable  method  of  triangulation  may  be  achieved  by  usin^ 
elementary  unitary  transformations  (QR-algorithm).   The  QR-algorithm  is 
defined  by  the  relations, 


Ak  =  \   Rk 

,  (3-6) 

Ak+1  =  \     Ak  \   =  Rk  \ 

where     R  is  an  upper-triangular  matrix, 
Q  is  a  unitary  matrix,  i.e., 
Q  =  Q_1 

The  unitary-triangular  decomposition  of  any  square  matrix  exists. 
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3.2  Theoretical  Background  (QR -Algorithm) 

(l)   For  any  matrix  A  (order  n),  there  exists  a  unitary  matrix 
Q  such  that 

\~  w 

where,  R,  is  an  upper  triangular  matrix  which  has  real  positive  diagonal 

elements,  and  Q,  =  IL  Nn  . .  •  W  (3-7) 

'  k    1  2      n  w  ' ' 

where 


N.  = 

l 


Ji-i  i   ° 

0    •   M. 


(3-8) 


M.  is  of  order  [n-(i-l)],  the  unitary  matrix  M  can  operate  on  any 
vector  b  such  that, 


M  b  =   |  |b  |  |  e 


(3-9) 


where 


|  To  |  |  is  the  Eulerian  norm  of  the  vector  b,  e,  =   {1,0,0,...0} 


M  is  an  elementary  unitary  matrix,  that  is  a  matrix  which 

differs  from  an  identity  matrix  at  most  in  one  principal  2x2  submatrix. 

Therefore,  if  the  vector  b  is  of  order  m,  M  =  T  T  -,  T„  Tn  (3-10) 

7  '  m  m-1       d     1 

where, 

t-n   0 0  t. 

11  lr 


T   = 
r 


1 0  0 

i  N   i 

i   \  I  i 

0 -1  0 


t  n   0 0  t 

rl  rr 
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The  elements  of  T  can  be  easily  expressed  in  terms  of  the  components 


r 


of  vector  b, 


<or-(t„  )r  =  [i .  JMi.2]  I 


ll7  rr'  L   |b_ 


p<r     » 


(3.11) 


(tn  r  =  -(t  j1  =  b  /  zib 

\    ]_r/  v    rl'  r'       '    p 

p<r 


Hence,  if  we  consider  that  the  matrix  A,  consists  of  the 

K. 


column  vectors  C.,  j  =1,  • • • ,  n,  therefore, 
J 


N1!^  N*  N^  A.  =  R.  (3-12) 

n   n-1       2   1   Tc    Is.  \~>        / 


where  N.  is  defined  by  Eq.  (3*8),  and  r..  =  |b.||,  where 

b.  =  {a..,  a.  .  . ,  .....  a   .},  i  =  1.  ....  n 

i   l  n'   i+l,i'     '   n,iJ'  ' 

(2)   Before  starting  this  step,  it  is  of  importance  to  prove 
that  Q  is  unique  if  A  is  non-singular. 


Proof:    Suppose  that 


A  =  Ql  Rl  =  Q2  R2 


therefore,  if  A  is  non-singular,  R  and  Rp  are  non-singular 

Rl  R~2  =   Ql"L  Q2  =  Ql  Q2 

therefore,  (R,  R?  )  is  unitary  also,  hence,  (R,  Rp  )  '"  =  (R-,  Rp  )    (3 '13) 

since  R  R„   is  an  upper-triangular  matrix,  Eq.  (3'13)  is  satisfied  if 

-1  -1      . 

and  only  if  R  R   is  a  diagonal  matrix,  therefore,  R  R?  =  I,  i.e. 

R1  =  R2,  hence  Q  =  Q   . 
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Step  (l)  is  completed  by  post  multiplying  (R  )  from  Eq.  (3-12) 

k 

by  Q  =  Nx  N2  ...  Nn  to  get  A    =  Rk  Q,  and  the  process  is  repeated  all 
over  again  until  all  the  elements  a. . ,  i  >  j  of  the  matrix  A  vanish  if 
the  eigenvalues  of  A  are  distinct.   Similar  to  the  LR-algorithm  the 
process  may  be  looked  upon  as  follows, 

if  Pk=QlQ2...Qk 


therefore, 


Sk  "  Rk  Rk-1  • • •  R2  Rl 


PkSk  =  QlQ2  ...  Q^  (Q^)  R^  ...  R2Rl 

=  Al  Ql  92    ••■    \J\-1  \-l>    \-2    •••    R2  Rl 


and  so  forth,    finally 

Tk  ~k  "   XJil 


P.   Sv   =   (Ajk  (3.1U) 


Since  P,  is  unitary  and  S   is  an  upper -triangular,  therefore  the  unitary- 
triangular  decomposition  of  a  non-singular  matrix  is  unique,  according 
to  the  lemma  proved  at  the  beginning  of  step  (2).   Francis  (1961)  proved 
the  convergence  of  P  S   =  (A)   and  showed  that,  "If  any  non-singular 
matrix  A  has  eigenvalues  of  distinct  modulus,  then  under  the  QR  trans- 
formation the  elements  below  the  principal  diagonal  tend  to  zero,  the 
moduli  of  those  above  it  tend  to  fixed  values,  and  the  elements  on  the 
principal  diagonal  itself  tend  to  the  eigenvalues." 

However,  for  eigenvalues  of  the  equal  modulus,  it  can  be  shown 
that  the  matrix  A,  becomes  split  into  independent  principal  submatrices 
coupled  only  in  so  far  as  the  eigenvectors  are  concerned.  Usually,  each 
of  these  submatrices  is  associated  with  groups  of  eigenvalues  of  the  same 
modulus. 
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(3)  The  amount  of  computations  that  the  method  involves  in 

3 
one  iteration  is  proportional  to  n  .  However  if  the  matrix  is  first 

reduced  to  the  upper  Hessenberg  form  by  any  stable  unitary  transformation 

technique,  i.e.,  Householder's  method,  the  amount  of  computations  become 

2 
proportional  to  n  .   Reducing  the  matrix  to  the  upper  Hessenberg  form 

has  some  advantages, 

(i)   If  any  element  a.  .  of  the  matrix  A  is  zero,  the 
matrix  can  be  partitioned  at  this  element  and  the 
rest  of  the  submatrices  can  be  operated  on  separ- 
ately so  far  as  the  eigenvalues  are  concerned. 

(ii)   The  upper-Hessenberg  form  is  preserved  under  the 
QR-transformation. 

(iii)   If  an  almost  triangular  matrix  A  is  non-singular 

and  has  non-zero  elements,  then  under  the  QR-trans- 
formation the  principal  diagonal  elements  of  A 
converge  to  the  eigenvalues  in  order  of  size,  if 
they  are  not  equal. 

(iv)  For  distinct  eigenvalues,  the  elements  (a   )   i  >  J, 

ij   k 

K     k 
below  the  principal  diagonal  converge  to   zero  as   , i,.    .      This   is  the 

^  \.} 
3 

same  rate  of  convergence  as  the  LR-transformation  which  is  rather  slow. 

Taking  a  closer  look  at  the  process  where,  for  example,  (a  )  approaches 

the  exact  value  of  the  eigenvalue  \    ,  it  is  found  that 

t>         n' 

ek+l  =  €k   *  ^ 

Vi 

where  ek  =  (a^  -  ^  and  \       <   1 

Vl 


To  speed  up  convergence,  the  origin  of  all  eigenvalues  is 

r   the  value  t     before 
3k 

after  the  iteration,  therefore 


shifted  by  the  value  £  before  an  iteration  and  then  shifted  back  again 
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n  -  bk 

€k+l   €k  *  X  -  t 

n-1   sk 


which  is  quite  an  improved  rate  of  convergence  for  t, ,  =  \   . 

k    n 

Therefore,  the  generalized  QR-transformation  becomes, 


Vi  =  \  \\ 


where  q*  (A)£  -  E  l)_  =  Ek 


In  Section  3  methods  of  choosing  the  origin  shifts  (\  are 

k 

discussed  for  "both  distinct  and  equal  modulus  eigenvalues  especially 
conjugate  complex  pairs.   Once  K     is  found  by  sufficient  accuracy,  the 
order  of  the  matrix  may  be  reduced  by  omitting  the  last  row  and  column, 


3.3  Numerical  Solution  for  QR  Algorithm 

The  numerical  solution  for  Q-R  algorithm  can  be  divided  into 
two  major  steps,  viz., 

(a)  Applying  Householder  's  technique  to  reduce  a  matrix 
to  an  upper  Hessenberg  form. 

(b)  Using  unitary  transformation  with  the  help  of  origin- 
shift  to  iterate,  i.e., 

A«  .  A 

A(K+1)  .  q(K)t  A(K)  Q(K) 

~(K)t  . (K)  produces  an  upper-diagonal  matrix, 
where  Q,     A 

Since  the  Householder's  technique  has  been  discussed 
in  detail  already,  therefore  we  will  only  concentrate  on  part  (b). 
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The  very  first  step  in  part  (b)  is  to  determine  the  value 
of  the  origin  shift.    To  find  this  value,  we  first  find  the  roots 
of  the  last  principal  2x2  submatrix  and  then  distinguish  the  follow- 
ing two  cases. 

(l)  Roots  are  real:   choose  one  that  differs  least  from 


the  last  diagonal  element  a  v  '    (before  kth  iteration) 

/ ,  \  nn  i \r\ 

and  call  this  root  X.        ,   then  compare  a.         with  the 
previous   one,   viz. a/        '(assuming  aA    '=0)    if 

1 


„(io  .  ^(k-D 


IW 


<2 


set  E^  =  K^\    otherwise  set  E^  =  0  where  E^ 
is  the  value  of  origin  shift  for  kth  iteration. 

(2)   Roots  are  complex:   we  do  one  of  two  things, 

(i)   Set  A.    =  |r|,  where  |r|  is  the  modulus 

of  the  complex  roots  and  then  compare  with 

(ii)   Do  double  origin-shift  by  setting  A.     =  R  and 

*  (k)    ^  -+v,  4.1,       •  \     (k-X) 

A.     =  R,  compare  with  the  previous  A. 

and  \p     ,  (assuming  a.    =  0,  a.   '  =  0; . 

After  the  value  of  origin  shift  has  been  found,  we  subtract 
all  diagonal  elements  from  this  value  and  then  start  iteration;  add 
the  value  back  after  each  iteration,  until  the  last  subdiagonal 
element  -*  0,    then  deflate  the  matrix  by  taking  out  the  last  row 
and  column,  the  last  diagonal  element  is  one  of  the  eigen-values. 

In  doing  QR  iteration,  we  simply  produce  a  series  of  2  x  2 
submatrices  to  operate  on  two  rows  at  a  time  on  matrix  A  and  eliminate  one 
subdiagonal  element  at  each  time  during  the  pre-multiplication. 
operate  on  two  columns  at  a  time  during  the  post-multiplication  accord- 
ing to  the  following  algorithm: 
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(i)      PRE-multiply: 


fi 


M-   V 


Nl 


1 
J 


r  r  r 

r 

r-   - 

- 

-   r 

r  r 

r 

r-   - 

- 

-  r 

P 

y 

y-  - 

- 

-  y 

a 

X 

x-   - 

- 

-    X 

a 

a-   - 

- 

-   a 

1 

a-   - 

V 

a 

-  a 

a  a 
a  a 

r  r 

r-   - 

- 

- 

r  r 

r  r 

r-   - 

- 

- 

r  r 

k 

y*y' 

- 

- 

yy 

x'x' 

- 

- 

x'x' 

a  a- 

- 

- 

a  a 

a- 

V 

- 

a  a 

\ 

a 

a  a 
a  a 

where  the  elements  r,  r'  and  k  are  those  of  the  final  triangular  matrix 
and  the  elements  a,  x,  and  0!  are  in  their  original  state.  The  elements 
in  the  two  rows  that  are  changed  are  given  by: 

k  =  (|p|2  +  \cc\2)2 
x ' =  ux  -  vy 
y'=  uy  -  vx 
u  =  p/k 
v  =  a/k 

(ii)   POST-multiply: 


a  x  y  r  r- 

a  x  y  r  r- 

a  x  y  r  r- 

k  r  r- 

r  r- 


N  _ 


xr 


f\ 


\1    -V 
V  [1 


\1 


a  a  x'y'r  r  - 

a  a  x'y'r  r  - 

a  x'y'r  r  - 

a   p  r  r  - 

r  r  - 

r  - 


-  r 

-  r 

-  r 

-  r 

-  r 

-  r 


Where  the  elements  r,  y  and  k  are  so  far  unaltered  from  the  triangular 
state,  and  the  elements  a,  x'  and  a  are  in  their  final  state.  The  two 
columns  that  are  changed  are  given  by: 
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a      =    vk 
P   =  jlk 

x'  =  \ix   +  vy 


yT   =  uy  -  vx 


The  higher  level  language  program  with  a.  detailed  flow 
chart  for  single  origin-shift  is  in  section  3«^+  and  3«5. 
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3«^  High-Level  Language  Program 


START: 


L00P1 : 


#BEGIN 

#LABEL 

#ARRAY 


#REAL 

#REAL 

#SET 
#READ 


#IF 


START,  LOOP1,  L00P2,  LOOP3,  LIP,  L2P,  L3P,  L^P, 
LOPM,  LOPE; 

#REAL  #FLOAT  #SHORT  #SKWD 

A[64,64],  P[64,64],  IM[6^,64],  UT[6U,l],  TEM[6U,l], 

V[61+],    \j[6k]; 

#FLOAT     #SHORT 

S2,    S,    T,    El,    TEMPI,    XI,   X2,    LI,   L2,    E2,    E,   K; 

#FIXED 

J,   H,   L,   Kl,   N,   M,    I; 

HE,   LL,   KK,    II,  IJ,    JI; 

A; 

J  -  2; 

HH: :    [J,    J+l,    ,    Gh] ; 

LL::    [l,    2,    ....,    J-l] ; 

KK:  :    [J+l,    J+2,    ,    Gk] ; 

S2  «-  #FOR   (H)   #SIM   (HH)   #DO     #SUM(A[H,  J-l]1'  2) ; 

S  «-  SQRT(S2); 

T  «-  S2   -  A[J,J-1]    x  S; 

#FOR   (L)      #SIM  (LL)      #DO     UT[L,1>0.; 

UT[J,1]   -  A[J,J-1]    -   S; 

#FOR   (Kl)  #SIM   (KK)     #DO     UT[K,l]   «-A[K,J-l]; 

P  «-  IM  -  TRANSPOSE (UT)   x  UT/T; 

A  <-  P  x  A; 

A  -  A  x  TRANSPOSE(P); 

J  «-  J+l; 

J  <  63  #THEN  #GOTO  LOOPl; 

N  *-  6h; 
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L00P2 


L00P3 


LIP 


L2P 


LOPM: 


II:  :    [1,    ?-,    ,   N] 

El  «-  0.; 

TEMPI  -  SQRT((A[N,N]   +  A[N-l,N-l] )t 2-k  x(A[N-l,N-l]   x 

A[N,N]    -  A[N,N-1]    x  A[N-1,N])); 
XI  «-    .5  x   (A[N,N]    +  A[N-1,N-1]    +  TEMPI); 
X2  -    .5   x   (A[N,N]    +  A[N-1,N-1]  -TEMPI); 
LI  «-  ABS(A[N,N]    -  XI); 
L2  «-  ABS(A[N,N]    -  X2); 
#IF  LI  <  L2     #THEN     #D0       #BEGIN 

E2  <-  LI; 
#GOTO  LIP; 
#END; 
E2  -  L2; 

#IF  ABS((E2-E1)/E2)    <  .5     #THEN     #D0     #BEGIN 

E  <-  E2; 
El  -  E; 
#GOTO     L2P; 
#END; 
E  *-  0.; 

#FOR      (I)      #SLM  (II)      #D0 
A[I,I]   -  A[I,I]    -   E; 
M  <-  1} 

JI: :    [M+l,M+2, ,M] ; 

K  ♦-  SQRT(A[M,M]    t2   +  A[M+1,M]  t2); 

U[M]   «-  A[M,M]/K; 

V[M]  <-  A[M+1,M]/K; 

#FOR   (I)  #SIM  (JI)   #D0  #BEGIN 

TEM[M,I]   <-  U[M]    x  A[M,I]    -   V[M]    x  A[M+1,I]; 

A[M+1,I]   -  U[M]    x  A[M+1,I]    +  V[M]   x  A[M,I]; 

A[M,I]   *-  TEM[M,I]; 

#MD; 

A[M,M]   «-  K; 

A[M+1,M]   *-  0.; 
#IF     M  >  N-l     #THEN     #GOTO     L3P; 
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M  *~  M+l; 
#GOTO     LOPM; 

L3P:  M  «-  1; 

LOPN:  IJ::    [1,2, ,M] ; 

A[M+1,M]  *-  V[M]    x  A[M+l,M+l]; 
A[ M+l, M+l]  <-  U[M]    x  A[M+l,M+l]; 
#FOR   (I)   #SIM   (IJ)   #DO  #BEGIN 

TEM[I,M]   *-  U[M]    x  A[I,M]    +   V[M]    x  A[l,M+l]j 

A[I,M+1]   -  U[M]    x  A[I,M+1]    +  V[M]    x  A[I,M]; 

A[I,M]  <-  TEM[l,M]j 

#EM); 

M  <-   M+l; 
#IF  M  >  N  #THEN  #GOTO  L^P; 

#GOTO  LOPN; 

L4P:  #FOR  (I)  #SIM  (II)  #DO 

A[I,I]  -  A[I,I]  +  E; 
#IF  ABS(A[N,N-1])  >  .OCOOOl  #THEN  #GOTO  L00P3; 
#PRINT  A[N,N]j 
N  <-  N-l; 
#IF  N  >  1  #THEN  #GOTO  L00P2; 
#PRINT  A[l,l]; 
#END 
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3-5  Flow  Chart 

TO  GET  UPPER  HESSENBERG  FORM 


H  «-  2 


2 
S   =  S2 


? 

2K~=  T 


m — 

S2-Z  [A(j,(H-l))]2 
J=H 


2K   _  s2  -  A(H,H-l)xS 


L—  x,2.f  •  •  •  jH-1 


UT(L)  «-  0 
UT(H)  «=  A(H,H-l)-S 


K=H+1,H+2,...,N 


UT(K)  -  A(K,H-1) 


I  is  64x6*4  

identity  matrix 


P^I -TRANS POSE ( UT ) xUT/ 2K' 


A  <-  PxA 
A  *-  AxTRANSPOSE(P) 


H  «-  H+l 


yes 


no 


© 
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FIND   ORIGIN-SHIFT  VALUE 


III)  • 


© 


N  «-  6k 


E1--0 


Xl4(a       +a     .        ,+SQRT((a       +a     ,        _  Y 
2X    n,n     n-l,n-l     ^     vv    n,n     n-l,n-l' 

-i+(a     .  *a       -a  *a     .      ))) 

n-l,n-l     n,n     n,n-l     n-ljrr'' 

X2«-i(a       +a     .        ,-SQRT((a       +a     _        .  )' 
2V    n,n     n-l,n-l  vv   n,n     n-1,11-1' 

-Ma     i        n*a       -a         n*a     .      ))) 
n-l,n-l     n,n     n,n-l     n-l,n 


E  -  E2 

El  «-  E 


LI 
L2 


a        -XI 

a       -X2| 
n,n        ' 


yes 


E2  «-  L2 


yes 
fc 


E  -  0 


A(l,l)-E-A(l,l) 


)   E2  <-  LI 


a       =A(N,N) 
n,n     v    '    ' 


1=1,2,..., N 
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PREMULTIPLY:   GET  UPPER  DIAGONAL  MATRIX 


M  *-  1 


K-SQRT(A2(m,m)+A2(m+l;m)) 


U(M) 
V(M) 


A(M;M)/K 

a(m+i,m)/k 


J=M+l,M+2, ,N 


a(m,j)«-u(m)*a(m,j)-v(m)*a(m+i,j) 

A(M+1,J)  <-U(M)*A(M+l,J)+V(M)*A(M,j) 


A(M,M)*^ 

a(m+i,m)*o 


Yes 


No 


*  M^M+1 


-  39 


POST -MULTI PLY  AND  TESTING 


© 


v 

M     -   1 

V 

a(m+i,mK-v(m)*a(m+i,m+i) 

a(m+i,m+i)«-u(m)*a(m+i,m+i) 

\ 

1 

a(j,mKu(m)*a(j,m)+v(m)*a(j,m+i) 
a(j,m+i)^u(m)*a(j,m+i)-v(m)*a(j,m) 

\ 

' 

M     «-  M+l 

i 

No 


M  >  N 


Yes 


A(l,l)-A(l,l)+E 


Yes 


N*-N-l 
PRT.A(N,N) 


J=l,2, ...,M 


1=1,2,..., N 


n  <  r 


No 


Yes 


PRT.A(1,1) 
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APPENDIX  A 

High-Level  Language  Program  For  Finding  MAX 

#BEGIN 

#LABEL   START,  LOOP,  LEND; 

#ARRAY   #REAL  #FLOAT  #SHORT  #SKWD 

A[6k,6k],   B[6k,    1],  Cl6k,    1]; 
#REAL  #FLOAT  #SHORT 

LIM,  CMAX; 
#REAL  #FIXED 

i; 

#SET  II,  JJ; 

START:         #READ  A 

II::  [1,2,...., 6k]; 

#FOR  (I)   #SIM  (II)   #D0 
B[I,1]  -  1.  ; 
LIM  ^  1.  ; 

LOOP:  C  -  A  x  B; 

CMAX  «-  C[l,l]; 

JJ::  [1,2,. ...,63]; 

#FOR    (I)      #SEQ   (JJ)      #DO     #BEGIN 

#LT     ABS(CMAX)    -  ABS(C[I+1, 1] )    >  0     #THEN     #GOTO     LEND 

#ELSE     CMAX  *-  C  [1+1,1]; 

LEND:  #END; 

#IF     ABS(ABS(CMAX)    -   ABS(LIM))    >   .00001     #THEN     #D0     #BEGIN 
LIM  «-  CMAX; 

#F0R    (I)      #SIM    (II)      #D0 
B[I,1]   -  C[I,1]/LIM; 
#G0T0     LOOP; 
#END     #ELSE 
#PRINT     CMAX; 


#END 
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READ  IN  A 


B[I,1]  -  I- 


< A  is  N  x  I  (6k   here) 


< 1=1,2,... ,M  (64  here) 


LIM  «-  1. 


C  <-  A  x 


FIND  MAX(C)  ->  CMAX 


yes 


B[I,1]  =  C[I,1]/LIM 


PRINT  CMAX 


STOP 


^3 


APPENDIX  B 


QP-Algorithm  With  Double  Origin  Shift 

When  a  matrix  has  a  complex-pairs  of  eigen-values,  the  double 
origin  shift  is  needed.  In  this  scheme,  we  combine  two  iterations  into 
one  with  a  pair  of  origin  shifts  since  single  iteration  may  not  be 

sufficient  enough  due  to  the  existence  of  a  complex -pair  eigen-values. 

th 
So  at  k   stage,  we  will  have 

A(k+2)  =  Q(k+i)*  Q(k)»  A(k)  Q(k)  Q(k+i) 

Let  W  =  Q(k)  Q(k+1)  and  P  =  (A(k)  -  E(k)l)  (A(k)-  E(k+l)l)  where 

E  '    and  E      are  the  pair  of  origin  shifts,  such  that  W*T  gives  an 
upper  triangular  matrix.  W*  here  is  unitary  and  W  =  N  N  ...  N 


N.  = 

i 


form: 


•i-1 


M. 

l- 


in  the  previous  discussion.   T  has  the  following 


x  x  x  x  x x 

x  x  x  x  x x 

x  x  x  x  x x 

x  x  x  x x 

x  x  x x 

X  X X 

N 

X 

x  XXX 


Let  the  first  column  be  r  then  N  should  be  defined  so  that 


N   r  = 
1   1 


rll l-el 
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Hence,  N   Av  J   N 


takes  the  form: 


X  X  X  X  X X  X 

X  X  X  X  X X  X 

X  X  X  X  X X  X 

X  X  X  X  X X  X 

X  X X  X 

X X  X 

V 

\ 

\ 

XXX 

X  X 


Where  the  new  matrix  has  few  more  elements  then  before,  and  this  is 
not  what  we  want  to  do,  therefore  the  general  technique  will  be  as 
follows: 


(l)   Perform  an  initial  transformation  (N   a 

( v\ 
A    with  double  origin  shifts. 


(k) 


N  )  on 


(2)  Reduce  the  resulting  matrix  to  almost  triangular  form 
by  a  method,  such  as  Householder's  to  obtain  A 


A  typical  stage  in  the  iteration  will  take  the  forms: 
Before  Householder's  transformation:   (ii)  After  Householder's  transformation: 


XXX  X'Y'Y'Y  Y  Y 

XXX  X'Y'Y'Y  Y  Y 

X  X  X'Y'Y'Y  Y  Y 

X  X'Y'Y'Y  Y  Y 

0  Z'Z'Z'Y  Y  Y 

0  Z'Z'Z'Y  Y  Y 

"  Z'Z'Z'T  T  T 

T  T  T 

T  T 


X 

X 

X 

Y 

Y 

Y 

Y 

Y 

Y 

X 

X 

X 

Y 

Y 

Y 

Y 

Y 

Y 

X 

X 

Y 

Y 

Y 

Y 

Y 

Y 

z 

Z 

Z 

Y 

Y 

Y 

Y 

z 

Z 

Z 

Y 

Y 

Y 

Y 

z 

z 

Z 

T 
T 

T 
T 

T 

T 
T 
T 
T 

T 
T 
T 
T 

where  3  rows  and  3  columns  being  affected.   X  and  T  are  those  of  the 

•  •  +  ■  1        *  «-   1    i.  i  «(k)     A(k+2) 

initial  and  final  matrices  A    and  A      res 

Y's  and  Z's  are  changed  by  the  transformation. 
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(k) 
Since  the  first  transformation  of  A    "by  N  is  similar  to 

the  subsequent  transformations,  we  can  set  up  the  initial  condition 

for  the  iteration  by  finding  r  ,  r_  ,  and  r   given  by 


rn  =  an  (an  "  a)  +  ai2  '  a2i  +  p 


r21  =  a21  (all  +  a22  "  a) 


r31  ==  a21  *  a32 


where 


a   =  E1  +■  E2 


E1  •  E2 


the  transformation  matrices  of  Householder's  method  are  of  the  form: 


N.  -  N.*  =  I  -  2t.t.* 

11  11 


t.  is  such  that   ti   =  1 
l  i  i   i  i 


t.  .  =  0  for  j  <  i  and  j  >  i+2 

1      Z10  i 


t.  . 


J20 


i,i+l    2kt 


li 


t.  . 


J30 


i,i+2    2kt, 


li 


2*i 


where  k  =(sign  Z1Q]  •  (Z1Q  +  Z2Q  +  Z^  )2 


where  Z 


rs 


Z10  Zll  Z12 

7  7  7 

20   21  22 

7  7  7 

30   31  32 


corresponds  to  the  Z's  of  the  figure  above, 
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In  order  to  reduce  the  number  of  multiplication  in  the  Householder's 

transformation  let  H.   =  [0, . . . ,0,1.x. ,Xo,0, . . . ,0]  =  t —  •  t. 

l  1   2  t..     l 

11 

obtaining  N.  =  I  -  (l  +  cp)  H.  H. 

i  11 

Z10         Z20  Z30 

*  =  T" '  Xl  =   k?Z^  '  X2  =  k^  and  a11  <  1 

Hence  at  i   stage,  the  operation  on  columns  a.,  a.  _  and  a.  „  is  [a., a.  ...a.  _] 

r -,   "I  i   i+l      i+2      i'  l+l'  H-2J 


[ai,  a1+1,  ai+2].[I  -  (l  +  cp) 


{a.  +  X  a    +  X2  a.+2] 


therefore,    a.  =  a.  -  n 
li 


Xl 
X2 


[1,  X1,  X2]]       Let  r)  =  (1  +■  cp)  . 


Vl  =  Vl  -  V11 


ai+2  =  ai+2  "X2^ 


Similarly  for  the  row  operation.   Z_0  and  Z~n  are  eliminated  and  replace 
Z  n  by  -  k  in  each  transformation. 

For  origin  shift,  at  k   stage,  we  need  to  consider  the  last 
2x2  sub-diagonal  matrix.   First,  find  the  two  roots  of  this  2x2  matrix, 

X        ,  X  which  differ  least  from  the  last  two  diagonal  elements 

(k-2) 

respectively.   Then  compare  with  the  previous  two  roots  X  and 

X  ,  and  calculate 

(k)     (k-2)        (k+1)     (k-1) 

I- ~r-^ 1  and  I- /  '■    v 1 

I    JJk) '     '    jjk+l) I 

if  both  >  |  then  set  E^  =  e^1'  =  0;  if  both  <  |  ,  then  set  E^k'  = 

.  (k)  „(k+l)   .  (k+l)   ,,    .     .  .  , ,  _(k)  ,  „(k+l)  ,  ,   ..     -, 
\v   ,  Ev    '    =  \s        ' ;    otherwise  set  both  Ev  andEv    '    to  be  the  real 

part  of  either  X^    '    or  X  whichever  corresponds  to  the  quantity 

<  p  .  Then  iterate  again  until  either  or  both  of  the  last  two  sub- 
diagonal  elements   are   near  zero,  then  deflate  the  matrix  accordingly; 
record  the  eigen-value  or  eigen-values;  find  new  origin  shifts;  iterate 
again  in  the  same  way  until  all  eigen-values  are  found. 
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''TRANQUIL  FOR  Q-R  ALGORITHM  WITH  DOUBLE  ORIGIN   SHIFT 

#BEGIN 

#LABEL    START,  LOOPl,  L00P2,  LOOP3,  ROWOP,  CLMOP; 

#ARRAY  #REAL  #FLOAT  #SHORT  #SKWD 

k[6k,6k],  v[6k,6k],  m[6k,6k],   ur[6i+,i]; 


START: 


LOOP1: 


L00P2: 


#REAL  #FLOAT  #SHORT 

S2,  S,  T,  El,  E2,  TEMPI,  XI,  X2,  LI,  L2,  EE1,  EE2, 
ESUM,  EPDT,  Rl,  R2,  R3,  K,  ALP,  PI,  P2,  AN,  RT1,  RT2; 
#REAL  #FIXED 

J,  H,  Kl,  L,  N,  I; 
#SET      HH,  LL,  KK,  II,  JJ; 
#READ      A; 

J-  2; 
HH::[J,J+l,...,6l4]; 
LL::[l,2,..,,J-l]; 
KK:: [J+l, J+2, . . . ,6k] ; 

S2  «-  #OR  (H)  #SIM  (HH)  #DO  #SUM(A[H,  J-l]t  2) ; 
S  -  SQRT(S2); 
T  <-  S2  -  A[J,J-l]*  S; 
#FOR  (L)  #SIM  (LL)  #DO  UT[L,1]  «-  0.; 
UT[J,1]  -  A[J,J-1]-  S; 

#FOR  (Kl)  #SIM  (KK)  #D0  UT[Kl,l]  *-  A[Kl,J-l]; 
P  «-  IM-TRANSPOSE  (UT)  x  UT/T; 
A  <-  PxA; 

A  ^  A  x  TRANSPOSE  (P); 
J  <-   J+l; 
#IF  J  <  63  #THEN  #GOTO  LOOPl; 
N  *-  6k; 

II::[1,2,...,N] 
El  <-  E2  <-  0.; 
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L00P3:  TEMPI  «-  SQRT((A[N,N]  +  A[N-l,N-l])t  2-l+*(A[N-l,N-l]* 

A[N,N]  -  A[N,N-l]*  A[N-1,N])); 
XI  *■  .5*  (A[N,N]  +  A[N-l,N-l]  +  TEMPI); 
X2  <-  .5*  (A[N,N]  +  A[N-l,N-l]  •-  TEMPI); 
LI  «-  ABS  (A[N,N]  -  XI); 
L2  «-  ABS  (A[N,N]  -  X2); 
#LE  LI  <  L2  #THEN  #D0  #BEGIN 

EE1  «-  XI; 

EE2  *-   X2;   #END  #ELSE  #D0  #BEGIN 

EE1  «-  X2; 

EE2  -  XI;      #END; 
#IF  ABS((EE1-E1)/EE1)   <   .5     #THEN 

El  «-  EE1     #ELSE 

El  <-  0.; 
#IF  ABS((EE2-E2)/EE2)   <   .5     #THEN 

R0WOP:  #FOR   (j)      #SIM(jj)     #D0     #BEGIN 

#LF   I  <  N   -   2     #THEN 

AN  <-  ALP*    (A[I,J]    f   PI*  A[I+1,J]    +   P2*  A[I+2,J]) 

#ELSE 

AN  <-  ALP*    (A[I,J]    +  PI*  A[I+1,J]); 

A[I,J]   <-  A[I,J]    -  AN; 

A[I+1,J]   «-  A[I+1,J]    -   PI*  AN; 
#IF  I  <  N-2     #THEN 

A[I+2,J]  <-  A[I+2,J]    -   P2*  AN; 
CLMOP:  #IF  I  <  N-2     #TKEN 

AN  -  ALP*    (A [J, I]    +■   PI*  A[J,I+1]    +   P2*  A[J,I+2]) 

#ELSE 

AN  +-  ALP*    (A[J,I]   +  PI*  A[J,I+1]); 

A[J,I]   <-  A[J,I]    -  AN; 

A[J,I+l]  «-  A[J,I+1]  -  PI*  AN; 
#IF  I  <  N-2  #THEN 

A[J,I+2]  *-  A[J,I+2]  -  P2*  AN; 
#END; 
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#BEGIN 


#IF  I  <  N-3  #THEN  #D0  #BEGIN 

AN  «-  ALP  *  P2*  A[l+3,I+2]j 

A[I+3,l]  -  -  AN; 

A[I+3,I+1]«--P1*  AN; 

A[I+3,I+2]  -  A[l+37I+2]  -  P2*  AN; 

#END; 

E2  <-   EE2  #ELSE 

E2  -  0.; 

ESUM  <-   E1+E2; 

EPDT  <-  E1*E2; 
#0R  (I)   #3EQ  (II)   #D0 

#LE   1=1  #THEN  #D0  #BEGIN 

Rl  «-  A[l,l]*  A(A[l,l]  -  ESUM)  +  A[l,2]  *  A[2,l]  +  EPDT; 

R2  «-  A[2,l]*  (A[2,2]  -  ESUM  +  A[l,l]); 

R3  -  A[2,l]*  A[3,2]; 

A[3,l]  *"  0;   #END  #ELSE  #D0  #BEGIN 

Rl  =  A[I,I-l]; 

R2  =  A[I+1,I-1]; 
#IF  I=N-1  #THEN  R3  -  0  #ELSE  R3  <-  A[I+2,I-l]; 

#END; 
#IF  Rl  >  0  #THEN 

K  <-  SQRT  (Rlt  2  +  R2t  2  +  R3t  2  )   #ELSE 

K  <-   -  SQRT  (Rlt  2  +  R2t  2  +  R3  t  2); 
#IF   K  ^  0  #THEN  #D0  #BEGIN 

ALP  <-   Rl/K  +  1.; 

PI  <-  R2/(RH-k); 

P2  -  R3/(Rl+k);   #END  #ELSE  #D0  #BEGIN 

ALP  <-  2; 

PI  <-  P2  <-  0.;  #END; 

JJ::  [1,1+1, ,N]; 

#IF  ABS(A[N,N-1])  >  .000001  #THEN  #D0  #BEGIN 

#IF  ABS(A[N-l,N-2])  >  .000001  #THEN  #G0T0  #L00P3; 

TEMPI  «_  SqRT((A[N,N]  +  A[N-l,N-l])t2-U*  (A[N--l,N-l]* 
A[N,N]  -  A[N,N-l]*  A[N-1,N])); 
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#END; 


#END 


RT1  *-   .5*  (A[N,N]  +  A[N-1,N-1]  +  TEMPI); 

RT2  «-  .5*  (A[N,K]  +  A[N-1,N-1]  -  TEMPI); 
#PRINT  RT1,  RT2; 

N  -  N-2;  #END  #ELSE  #D0  #BEGIN 
#IF  ABS(A[W-l,N-2])  <  .000001  #THEN  #D0  #BEGIN 

#PRINT  A[N,N],  A[N-l,N-l]; 

N  <-  N-2;   #END  #ELSE  #D0  #BEGIN 

#PRINT  A[N,N]; 

N  «-  N-l;  #END; 
#END; 

#IF  N  <  1  #THEN  #G0T0  L00P2; 
#ST0P; 
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TO  GET  UPPER  HESSENBERG  FORM 


H  «-  2 


2 
S  =   S2 


2 

2K   5  T 


2 


"FT 

I  [A(J,(H-1).)]! 

J=H 


2K2  *-  S2  -  A(H,H-l)xS 


L=1,2,...,H-1 * 


K=H+l,H+2,...,N-> 


UT(L)  *-  0 
UT(H)  <-  A(H,H-1)-S 


UT(K)  <-   A(K,H-1) 


I  is  6kx6k 
Identity 


P  <-  I -TRANSPOSE  (UT)xUT/2K' 


A  -  PxA 
A  «-  AxTRANSPOSE(P) 


H  <-  H  +  1 


yes 


52  - 


FIND  ORIGIN  SHIFT 


d> 


El  «-  0  « 


E2  «-  0 


N  «-  6U 


El  <-  E2  «-  0 


Find  two  roots  of  the 

last  2x2  matrix  call 

XI,  X2 


LI 
L2 


a   -XI 
n,n 

a  '  -X2 


ESUM  «-  E1+E2 
EPDT  «-  E1*E2 


EE1  -  XI 

EE2  «-  X2 


El  «-  EE1 


E2  «-  EE2 
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SET  UP  AND  DO  ROW  AND  COLUMN  OPERATIONS 


© 


i=2,3,-..,N-2-» 


i=N-l 


Rl=a.  .  . 


R2=a, 


Rl=a    (a   -ESUM)+a   *a   +EPI 
XfX      1,1         ±,d      ^->1 


!,l(a2r 


-ESUM+a   ) 


R>a2,l  *  a3,2 


I   2   2   2 

K  =  -jRl  +R2  +R3 


yes 


ALP  «-  Rl/K  +  1 
PI  -  R2/(R1+K) 
P2  «-  R3/(R1+K) 


ROW  OPERATION 


COLUMN  OPERATION 


2   2   2 
=  jRl  +R2  +R3 


ALP  4-  2 
PI  <-  P2  ♦-  0 
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END  TEST 


PRT.A(N,N) 
N  *-  N  -  1 


no 


III 


no 


yes 


STOP 


Find  two  roots  of  the 

last  2x2  matrix  call 

RT1,  RT2 


PRT.  RT1,  RT2 
N  <-  N  -  2 
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